///////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//-------------------------------- Ebrahim Foulaadvand, 14 May 2012 ---------------------------------------
//
// --- The routine "ChemicalMaster" evaluates the temporal evolution of the average number of a
// --- species in a chemical reaction.
//
//     tau=time step, m=initial number of type X molecules, k=reaction rate
//
//     N[t]=number of X type molecules at timestep t.
//
///////////////////////////////////////////////////////////////////////////////////////////////////////////


#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
#include<iomanip.h>


using namespace std;


double  tau=0.0000001,k=1,Avnum;

int M=151,m=150,T=2000,t;

main()

{

ofstream  fileAvnum ("Avn-t k=1 m=150.plt");

vector<double> P(M+2,0),Avn(T+1,0);


// ---------------------------------- initial conditions --------------------------------------------

for(int n=0;n<=M;n++){ P[n]=0; }

P[m]=1;

Avn[0]=m;

//---------------------------------------- Eule Method-----------------------------------------------

for(int k=1;k<=T;k++){

Avnum=0;

P[0]=P[0] + tau*k*P[2]; //cout<<"P[0]= "<<P[0]<<endl;//getch();


for(int n=1;n<M-1;n++){

P[n]=P[n] + tau*k*0.5*((n+2)*(n+1)*P[n+2]  - n*(n-1)*P[n]); //cout<<"P["<<n<<"]= "<<P[n]<<endl;getch();

Avnum+=n*P[n];

}

P[M-1]=P[M-1] + tau*k*0.5*( (M+1)*M*P[M] - (M-1)*(M-2)*P[M-1] );

Avnum+=(M-1)*P[M-1];


P[M]=P[M] + tau*k*0.5*( (M+2)*(M+1)*P[M] -  M*(M-1)*P[M] );

Avnum+=M*P[M];

fileAvnum<<k*tau<<" "<<Avnum<<endl;

//cout<<"Avnum= "<<Avnum<<endl;//getch();

//cout<<"P["<<m<<"]= "<<P[m]<<endl;getch();


} // end k

cout<<"end."<<endl;getch();

return 0;

}